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ABSTRACT 

The power spectrum P{k) oc /c"’ with n = —2 is close to the shape of the 
measured galaxy spectrum on small scales. Unfortunately this spectrum has 
proven rather difficult to simulate. Further, 2-dimensional simulations have 
suggested a breakdown of self-similar scaling for spectra with n < —1 due to 
divergent contributions from the coupling of long wave modes. This paper is 
the second (numerical) part of our investigation into nonlinear gravitational 
clustering of scale-free spectra, in particular to test the scaling of the n = —2 
spectrum. 

Using high-resolution N-body simulations we hnd that the n = —2 power 
spectrum displays self-similar scaling. The phase shift of Fourier modes of the 
density show a dual scaling, self-similar scaling at early times and a scaling 
driven by the bulk velocity at late times. The second scaling was shown 
analytically to be a kinematical effect which does not affect the growth of 
clustering. Thus our analytical and N-body results verify that self-similarity 
in gravitational clustering holds for —3 < n < 1. The N-body spectrum 
is also compared with analytic htting formulae, which are found to slightly 
underestimate the power in the nonlinear regime. The asymptotic shape of 
the spectrum at high-fc is a power law with the same slope as predicted by the 
stable clustering hypothesis. 


Subject headings: cosmology: theory — large-scale structure of universe — 
galaxies: clustering — galaxies: formation 


1. Introduction 
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The self-similar scaling of density perturbations with scale free initial conditions in 
a spatially flat universe is a useful theoretical tool in studying structure formation. It 
has been widely used to study gravitational clustering in cosmology and has been tested 
by several studies using N-body simulations. However self-similar scaling for scale free 
initial spectra P\in{k) oc k"" has not been adequately demonstrated for n < —1 because the 
requirements of dynamic range get increasingly difficult to meet as n gets smaller. Indeed 
results of some two dimensional studies suggest a breaking of self-similar scaling for n = —2 
in three dimensions. Analytical analyses have been limited to the observation that the 
linear peculiar velocity held diverges for n < — 1, but the linear density contrast does not 
diverge provided n > —3. This would suggest that while there may be formal problems 
with establishing self-similarity for n < —1, in practice it should hold as long as n > —3. 

In an earlier paper (Jain & Bertschinger 1996 - Paper I), we have analyzed the 
dynamics of the coupling of long wave modes by analytical techniques to address the 
question of whether self-similar scaling is broken for n < — 1. On the basis of the nonlinear 
growth of density perturbations we concluded that self-similar scaling does not break 
down provided n > —3, consistent with the linear theory expectation. We also identihed 
statistical measures which would scale differently from self-similar scaling owing to the 
kinematic effect of large-scale bulk flows. These measures are related to the phase of the 
Fourier space density held, but do not affect the growth of its amplitude. 

In this paper we examine the scaling behavior of the n = —2 spectrum using high 
resolution three dimensional N-body simulations. Testing self-similarity using simulations is 
crucial because the analytical demonstration of a scale-free dynamics provides no guarantee 
that the similarity solution will actually hold, starting from some general initial conditions. 
Moreover, the analytical results of Paper I relied on an assumption about the initial helds, 
as well as on using the huid limit approximation. The N-body simulations do not rely 
on these approximations, though as we shall see their hnite resolution introduces other 
departures from exact self-similarity. Hence by combining our analytical and numerical 
results, we attempt to provide a complete analysis of the scaling behavior of scale free 
spectra for —3 < n < 1. 

We shall consider the similarity properties of gravitational dynamics in a zero-pressure 
Einstein-de Sitter cosmology. An Einstein-de Sitter universe refers to the model with the 
cosmological density parameter H = klmatter = 1 and zero cosmological constant, so that the 
universe is spatially flat. The gravitational interaction also does not pick a special length 
scale. Further let the initial power spectrum be a power law, Piin{k) oc k'^, over length 
scales of interest. In such a case there is no preferred length scale in the universe, and the 
evolution of structure is expected to be self-similar in time. 
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The resulting self-similar scaling of characteristic length scales a:, and wavenumber 
scales k is (Peebles 1980, Section 73): 

Xss{t) oc ; kss{t) ~ Xss{t)~^ oc (1) 

where a{t) is the expansion scale factor. The above scaling is most simply derived by 
requiring consistency with the growth of the rms linear density contrast, smoothed on a 
scale Xss = k~g , and related to the power spectrum as: 



where W{kx) is the smoothing window function. Requiring that the scale at which 
(Sp/p) = 1 scales self-similarly gives the result in equation (|I]). 

Early studies of self-similar evolution in cosmology include those of Peebles (1974); 
Press & Schechter (1974); Davis & Peebles (1977); and Efstathiou & Eastwood (1981). 
Efstathiou et al. (1988) tested self-similar scaling in N-body simulations of scale free spectra 
with n = —2, —1,0,1. They examined the scaling of the correlation function ^{x,t), and 
of the multiplicity function describing the distribution of bound objects. They verihed the 
predicted scaling for both statistics, and found consistency with the picture of hierarchical 
formation of nonlinear structure on increasingly large length scales. Their results for 
n = —2 did not match with the self-similar scaling as well as the other cases. Bertschinger 
& Gelb (1991) used better resolution simulations to address these questions and also found 
similar results. These authors concluded that the reason for the weakness of the n = —2 
results was the hnite size of their simulation box, as the n = —2 case has more power on 
large scales and therefore requires a larger box-size to approximate the inhnite volume limit 
with the same accuracy as larger values of n. 

More recently Lacey & Cole (1994) have examined the self-similar scaling of the 
number density of nonlinear clumps for scale-free spectra. Their results indicate that 
self-similar scaling works reasonably well for the statistics they measure, even for n = —2. 
Jain, Mo & White (1995) and Padmanabhan et al. (1995) reach the same conclusion for 
the correlation function and power spectrum. Colombi, Bouchet & Hernquist (1995) have 
verihed self-similarity for higher moments of the density. However, in all the above studies 
the self-similar scaling of a particular statistic is tested without making a comparison with 
alternate scalings. As demonstrated in Section 4, this makes it difficult to distinguish the 
effect of limited numerical resolution from a real breakdown of self-similar scaling. 

The N-body results of Ryden & Gramann (1991), and Gramann (1992) suggest that 
the n = —2 case is different from n > — 1 for a fundamental reason. They studied n = — 1 
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simulations in two dimensions, which are the analog of n = —2 in three dimensions, 
and examined the scaling of the phase (Ryden & Gramann 1991), and then both phase 
and amplitude (Gramann 1992) of the Fourier transform of the density field. The 
scaling was found to be different from the standard self-similar scaling. Gharacteristic 
wavenumber scales, instead of following the self-similar scaling, given in two dimensions 
by kss{t) oc oc showed the scaling k oc Other studies in two 

dimensions also suggest that a transition in nonlinear evolution occurs at n = —1 (Klypin & 
Melott 1992, and references therein). Motivated by Gramann’s results we had re-examined 
the n = —2 simulation presented in Bertschinger & Gelb (1991) and found that the results 
were ambiguous, and that a bigger simulation would be needed to provide a definitive 
answer. 

In this paper we use the n = —2 case as representative of the range —3 < n < — 1 and 
use high resolution N-body simulations to test for self-similar scaling. Section 2 provides the 
analytical background and motivates the need to test the scaling of spectra with n < — 1. 
The N-body simulations used to test self-similar scaling are described in Section 3, and 
effects of finite numerical resolution are discussed. In Section 4 results for the scaling of the 
power spectrum are presented and self-similar scaling is compared with an alternate scaling. 
In Section 5 a different statistic is used which measures the scaling of the amplitude as well 
as the phase of the Fourier space density. We discuss subtleties in measuring the evolution 
of the phase, and provide a kinematical interpretation of the results for the phase shift. We 
conclude in Section 6 by discussing the results of this paper and Paper I. 


2. Analytical Background 


This section summarizes the formalism and motivation presented in Section 2 of 
Paper I. We use comoving coordinates x and conformal time dr = dt/a{t) to write the 
nonrelativistic cosmological fluid equations as 
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+ V- [(1 + (5)F] = 0 , 

(3a) 

fF-vW = —-F—V0 , 

(3b) 

^ ^ a 

= AnGa^pd , 

(3c) 


where d = da/dr, 6 {x,t) = {p{x,t) — p{t))/p{t), (f>{x,T) is the perturbed gravitational 
potential, and v{x,t) = dx/dr is the proper peculiar velocity. We assume an Einstein-de 
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Sitter (f2 = 1) universe, with a oc oc r^. We will also assume that the initial (linear) 
density fluctuation fleld is a Gaussian random fleld. 

The fluid equations (|^) are strictly valid only on scales large compared to the nonlinear 
clustering scale. On scales much smaller than the nonlinear scale, the intersection of 
particle trajectories leads to a complicated, anisotropic stress tensor and thus invalidates 
a fluid description which requires the density and velocity to be smooth, single-valued 
functions of position. Whereas the analytical analysis of Paper I used the zero-pressure fluid 
approximation and the continuum limit, the N-body simulations do not as they integrate 
the equations of motion for the particle trajectories. 

To quantify the amplitude of fluctuations on various scales it is preferable to work with 
the Fourier transform of the density fluctuation fleld, which we define as 

Hk,T)= I-^e-'‘-S(x,T) . ( 4 ) 

The power spectrum (power spectral density) of S{x,t) is defined by the ensemble average 
two-point fnnction, 

{5{ki, T)S{k 2 , r)) = P{ki, t) dnih + h) , (5) 

where (5 d is the Dirac delta fnnction, required for a spatially homogeneous random density 
fleld. For a homogeneous and isotropic random fleld the power spectrnm depends only on 
the magnitnde of the wavevector. The contribution to the variance of S{x, r) from waves in 
the wavevector volnme element d^k is P{k, T)d?k. The antocorrelation fnnction is defined 
as ^{x,t) = { 6 {xi,t) 6 {x 2 ,t)), where \xi — X 2 I = x. It can be easily verified that ^(x,r) is 
the Fonrier transform of P{k,T). 

The analytical approach developed in Paper I can now be snmmarized as follows. For 
a scale free initial spectrnm P\in{k) oc fc", the rms density contrast in linear theory is given 
by equation (^. As fc —> 0 the window fnnction generically (e.g. the top-hat or Gaussian 
window fnnctions) has the limiting form W{kx) 1. Hence for n < —3 the linear rms 
density contrast on any scale receives divergent contribntions from the fc —>■ 0 part of the 
spectrnm. Analogous to the rms density contrast one can define the bnlk flow velocity, 
which is the average, smoothed peculiar velocity on scale x. It is obtained by using the 
peculiar velocity power spectrum Pun^ in equation (^ instead of Pun. From eqnation (pa|) 
it follows that 5 = —V • v in linear theory, or P\\-av{k) = a?/a? P\\-a{k)/k!^ . Thus the bulk 
velocity is given by 

vI{x,t) = 0^4:71 J dkPiia{k)W‘^{kx). ( 6 ) 

For n < —1, Vb{x,T) diverges dne to the fc —>■ 0 contribntion. Since the nonlinear flnid 
eqnations conple the density and velocity, this divergence manifests itself in the nonlinear 
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terms in the fluid equations. 

If this divergence were to have a dynamical influence in the growth of the density 
contrast it would lead to a breakdown of self-similar scaling, because it introduces a new 
scale in the system. Setting Vb{x) = constant < c, with a finite lower limit in the integral 
in equation (H) gives the scaling a: oc a for all n < — 1. The standard self-similar scaling 
of equation (|l]) then holds only if n > —1. This, coupled with the numerical results 
of Gramann (1992) described in Section 1, provided the motivation for our analytical 
examination of spectra with n < — 1 in Paper I. The results of Paper I are summarized in 
Section 5. In the following sections we will analyze the scaling of the n = —2 spectrum in 
N-body simulations to test whether it follows the velocity scaling a: oc a or the self-similar 
scaling of equation (|^): x (x. a?. 


3. The N-body Simulations 

N-body simulations provide a powerful means for testing the self-similar scaling of scale 
free spectra. The deeply nonlinear regime is accessible in these simulations, thus offering 
the possibility of measuring the complete similarity solution. N-body simulations have 
limited dynamic range, but they do not rely on any approximations of the kind made in 
Paper I. Since the equations of motion for individual particles are integrated, the fluid limit 
approximation is not required either. Therefore they provide a complementary technique to 
the analytic approaches of Paper I. 

The simulations used for the self-similar scaling analysis are three dimensional 
particle-particle/particle-mesh (P^M) simulations. The larger of our n = —2 simulations 
had 256^ particles and a Plummer force softening parameter e = 1/5120 L, where L is 
the box-length. It was performed on the CM-5 supercomputer at the National Center for 
Supercomputing Applications using the parallel code of Ferrell & Bertschinger (1994, 1995). 
This simulation, labeled A in Table 1, is used for the main results of this paper. To test 
the effects of numerical resolution at high-/c we have compared the results with a smaller 
P^M simulation (labeled B), performed on a workstation, which used half as many particles 
and twice the force softening length. This simulation used a variable PM mesh as shown in 
table 1. We have also used a set of three PM simulations (labeled C), performed on the 
CM-5, to estimate the effects of statistical fluctuations at low-fc. Table 1 gives the relevant 
parameters of all the simulations used. 

The initial positions and velocities for the particles were generated using the Zel’dovich 
approximation, and assumed the linear power spectrum and Gaussian statistics. Outputs 
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of the particle positions were stored at time intervals corresponding to an increase of a by a 
factor of 2^/^. Energy conservation was better than 10“^ of the total gravitational energy 
at all times for simnlation A and was similar for the other simnlations. 

At each of these ontpnt times the particle positions were interpolated onto a 1024^ grid 
nsing the triangnlar-shaped-clond (TSC) scheme (see Gelb & Bertschinger 1994 for details) 
to get the real space density. The density held 6{x, a) was then fast Fonrier transformed to 
get S{k, a). S{k, a) is a complex nnmber at each fc, and is therefore represented by a real and 
imaginary part. The valnes of k are represented by a three dimensional vector {k^, ky, kz), 
with each of k^, ky and kz being integers ranging from —511 to 512. In the remainder of 
this section the magnitnde of k will be given in nnits of 211 jL, so that the modes with 
wavelength eqnal to L have k = 1. The power spectrnm was compnted as the mean of 
\ 6 {k, a)p in bins of k with width 6 k = 1. The power spectrnm was corrected for the effects 
of shot noise and of convolntion with the TSC interpolating window nsed to compnte the 
real space density. Fignre 1 shows the power spectrnm for simnlation A at hve different 
valnes of a, corresponding to valnes of the nonlinear wavennmber (dehned in eqnation (|^) 
below) kni{a) = 58, 29,14.5, 7.25. Also shown as the dotted line is the linear spectrnm at 
the same times. 

In an N-body simnlation there are many departnres from the idealization of an inhnite, 
continnous flnid. The limitations common to any N-body simulation are the hnite size of 
the box and the discrete nature of particles. Additional scales are introduced due to the PM 
mesh, force softening and the density interpolation grid. As discussed in Section 1 of Paper 
I, the presence of these scales is not a fundamental drawback; it only means that one must 
be careful to ensure that the range of scales used to study the scaling in time is sufficient 
to allow for intermediate asymptotic self-similarity to set in. In the following discussion we 
shall focus on the largest of our simulations, labeled A in Table 1. 

An important scale that we shall refer to is the nonlinear scale, at which the fluctuation 
amplitude is ~ 1. We dehne a nonlinear wavenumber as: 

fknl 

/ d^kPUa,k) = l. (7) 

Jo 


Table 1: Parameters of the N-body simulations 


Simulation 

knl/{211/1) 

^part 

Softening e/L 

PM Mesh 

Timesteps 

A: n = —2, P^M 

5250 - 7.25 

256^ 

1/5120 

512^ 

1331 

B: n = —2, P^M 

2060 - 4.02 

128=^ 

1/2560 

256^ - 432^ 

1007 

C:n= -2, PM 

5250 - 7.25 

256^ 

— 

512^ 

692 












For Piin(a, A;) = Aa^k~‘^, equation gives kni{a) = (47rAa^)“F Thus kni{a) becomes 
smaller at late times since the nonlinear length scale increases with time. For scale free 
spectra kni{a) provides a useful time variable, for example, in comparing simulations with 
different starting amplitudes. We now turn to the main numerical limitations in our 
simulations. 

(i) Initial amplitude. The initial amplitude used in our n = —2 simulations A and C 
is such that at a = 1 the dimensionless power on the Nyquist frequency of the particle 
grid, k = 128 is 47rk^P{k) = 0.024. This corresponds to having kni = 5250 at the initial 
time. The n = —2 spectrum needs to be started with very small initial amplitude in order 
for clustering on small scales to develop self-similarly. As pointed out by Lacey & Cole 
(1994) the standard prescription of Efstathiou et al. (1985) of having power on the Nyquist 
frequency of the particle grid equal to the white noise level is too large and produces 
signihcant departures from self-similar evolution. 

(ii) Force softening. Force softening is implemented in computing short range forces 
in the particle-particle part of the P^M simulations by using the Plummer force law: 

F{r) = Grnf j (r^ -|- e^). The parameter e is given in Table 1 for the various P^M simulations. 
The distance at which r‘^F/{Gm^) = 1/2 is ~ 1.3e. This distance is ~ L/4000 for simulation 
A: it is smaller by almost a factor of 4 than the main limiting factor on small scales, the 
spacing of the 1024^ grid used for interpolating the density (as shown below in point (iii), 
the grid effects are in fact significant up to twice the grid spacing). This is important, since 
the force softening affects the dynamical evolution and can therefore artihcially suppress the 
clustering on scales several times larger than e. In addition, the hnite number of particles 
and the hnite size of the timesteps used limit the numerical resolution on small scales. For 
the parameters in our simulations these affect the spectrum on about the same limiting 
scale as the force softening. Hence their effect is also far less important than that of the 
FFT grid, aside from the PM simulations whose force softening was relatively large. 

(iii) FFT grid. The density held was fast Fourier transformed to get 5{k). This involves 
interpolation on to a grid to get the real space density held, and thus imposes a minimum 
length scale below which power cannot be measured. Figure 2 shows a comparison of the 
power spectrum computed for simulation B using two diherent grid sizes: 512^ and 1024^. 
Both spectra have been corrected for the ehects of convolution with the TSC window. It 
shows that the power for k up to 1/2 the Nyquist frequency of the FFT grid is accurate 
to better than 20%. However this is only a consistency check, as a conservative estimate, 
we therefore use a maximum k = 256 for the 1024^ interpolation and FFT grid in our 
self-similarity analysis. 

(iv) Statistical fluctuations at low-k. At low wavenumbers the hnite number of modes 
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present in a given bin in k causes statistical fluctuations in the power spectrum. This 
is evident in Figure 1 which shows that dX k = 1 and 2 there are signihcant differences 
between P{k) and P\\n{k) even at early times. By checking the power spectra of the three 
different PM runs of simulation C we found that as the simulation evolved the power in 
realizations with lower initial power at A; = 1 and 2 was suppressed out to A; ~ 5 due to 
nonlinear mode coupling. For the latest time output used in our analysis, at A; > 5 these 
fluctuations are smaller than 20%. For the results shown in Section 4 we plot the power 
spectra from k = 3 — 256. 

(v) Range of kni{a). At late times the nonlinear wavenumber knfa) approaches the 
size of the box. This sets a maximum value of a beyond which the absence of modes with 
wavelength larger than the box-size inhibits the growth of the spectrum at late times, 
because the coupling of these modes would have enhanced the power at high-A; (Jain 

& Bertschinger 1994). To test this effect, simulation B was evolved to a maximum a 
corresponding to kni = 2. Using this simulation we hnd that for the n = —2 spectrum the 
power at A; % 10 is indeed underestimated for kni{a) % 4. At later times, for example with 
kni{o.) = 2.8, the power spectrum at k extending up to A; ~ 100 is signihcantly suppressed 
relative to the amplitude required to satisfy self-similarity. This problem is particularly 
severe for the n = —2 spectrum because of the large amount of power on low-A; relative to 
spectra with larger n, which in turns leads to a stronger nonlinear coupling with high-A; 
modes. It is the principal reason for the requirement of a large dynamic range to test the 
scaling properties of this spectrum. 

Our choice of kni{amax) = 7.25 for simulation A is adequate to get dynamically 
accurate power spectra at low-A;, and also provides a sufficient range in a to test self-similar 
evolution. At small a we are limited by the fact that the shot noise amplitude which we 
subtract from the measured power spectrum exceeds the measured spectrum at high-A; for 
a < 8 {kni = 82). This determines the earliest output used in our analysis - for simulation 
A this had kniiamin) = 58. 

(vi) Accuracy of PM simulations. A hnal point, regarding the adequacy of PM 
simulations to follow the scaling of the n = —2 spectrum is made in Figure 3, where the 
power spectrum of the P^M simulation A is compared with the PM simulation C. These 
simulations have identical number of particles, the same initial Fourier modes, and are 
started with the same kni{a). Both spectra are plotted upto 1/2 the Nyquist frequency 
of the FFT grid, and can therefore be legitimately compared over the full range of k for 
which they overlap. The figure shows that at late times the power in the PM simulation 
is suppressed at wavenumbers significantly smaller than the Nyquist frequency of the PM 
grid, k = 512. At the last output time shown, corresponding to kni{a) = 7.25, the power 
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at A; > 50 is suppressed by > 20%. The suppression increases to > 50% aX. k > 80, i.e. at 
k > l/6th the Nyquist frequency of the PM grid. Thus the commonly cited limit of the 
spatial resolution of PM simulations as being less than twice the grid spacing considerably 
overestimates the spatial resolution. Our results are at odds with previous tests of PM 
simulations in the literature (e.g. Splinter, Melott & Shandarin 1997). The difference in 
these tests is the statistical measures used as well as the degree of nonlinearity, since the 
discrepancy appears to occur only in the regime of very high over-densities (the last output 
time shown). 


4. Scaling of the Power Spectrum 


The self-similar solution for the power spectrum and autocorrelation function is 
(Peebles 1980, Section 73), 

P{k,r) = a^°‘kQ^P{ka'^/ko) ; ^(x, r) = ^{x/xoap , (8) 

where a = 2/(3 -t- n); kQ, xq are constants which must be determined from the initial 
conditions; and P, ^ are unspecihed dimensionless functions. It is easy to verify that the 
linear spectrum Piin(A:, r) oc a?k^ is consistent with the functional form of equation (^. 
This functional form provides a strong constraint on the evolution of the spectrum, as P 
is no longer an arbitrary function of two independent variables k and a, but only of the 
variable fca“. 

The results for the scaling of the power spectrum are shown in Figures 4 and 5. The 
self-similar scaling given by equation (H) with n = —2 is compared with the velocity scaling. 
The two scalings are distinguished by the way in which characteristic wavenumbers kc scale 
with time: kc oc a“", with a = 2 for self-similar scaling, and a = 1 for the velocity scaling, 
as discussed at the end of Section 2. In Figure 4 we show the measured power spectrum at 
eight different times corresponding to a range of kni{a) = 58 — 7.25 (in units of 2'7i/L). Also 
shown as the dotted curves is the power spectrum from the fourth from last time, re-scaled 
self-similarly in the left panel, and according to the velocity scaling in the right panel. 
The results show that the measured spectra are consistent with self-similar scaling and 
show departures from the velocity scaling which causes the spectrum to grow too slowly 
at high-fc. It should be noted that there are apparent departures from self-similar scaling 
at low-fc. We believe these arise due to the statistical fluctuations resulting from the small 
number of modes present in the low-A: bins. These can be seen even in hgure 1 where the 
power at the lowest k does not agree with the linear spectrum. At the latest times however 


-11 - 


this discrepancy extends to fc ~ 10. A possible explanation is the missing power from 
long-wave modes, which, as one might expect from the action of nonlinear mode coupling, 
appears to affect the spectrum at increasing k at later time outputs. However this issue 
merits further exploration as it contradicts our own estimate of fc > 4 being an adequate 
cut-off at the low-/c end. 

Figure 5 shows the same feature in a slightly different way: the power spectrum curves 
at the eight different times are re-scaled so that they would lie along a single curve if they 
followed the predicted scaling. Again, the left panel which uses self-similar scaling shows 
that the curves almost coincide. The right panel showing the velocity scaling has significant 
differences between the different curves even at intermediate k on which the spectrum is 
measured very accurately. In figures 4 and 5 the agreement with the self-similarly scaled 
spectra is typically within 20%, and better than 40% over nearly the entire range (excluding 
only the lowest k where fluctuations occur due to the small number of modes). For the 
velocity scaling, the agreement is typically at the 50% level and is about a factor of 2 if 
one uses the earliest and latest output times. The curve traced out by the scaled N-body 
spectra in the left-panel of hgure 5 give the empirical P of equation (|). 

To highlight the numerical difficulty of demonstrating self-similar scaling for the 
n = —2 spectrum, in Figure 6 we show the results from the smaller of our P^M simulations, 
labeled B in Table 1. The figure shows the scaling of the spectrum in exactly the same 
way as Figure 5. On comparing the two panels in Figure 6 it is evident that the dynamic 
range of simulation B is not adequate to discriminate between the self-similar and velocity 
scalings. The scatter in the different curves due to numerical effects at low and high-/c 
(which overlap with intermediate-^ spectra at other times because of the re-scaling made in 
the plots) is comparable to the difference in the self-similar and velocity scalings. It is also 
clear that without comparing the two scalings, it is very difficult to distinguish the effects 
of limited dynamic range from real departures from the underlying scaling. 

The resolution of simulation B shown in Figure 6 is comparable to the simulations used 
by Lacey & Cole (1994) and Jain, Mo & White (1995). Lacey & Cole used a P^M simulation 
with 128^ particles and e = L/4160, while the simulation used by Jain, Mo & White was 
simulation B itself. Colombi, Bouchet & Hernquist (1995) used a treecode simulation 
with 64^ particles and a spatial resolution somewhat larger than that of simulation B. 
The resolution limitations of the PM simulation C are even more severe than those of 
simulation B. At late times the power at high-fc is significantly smaller than that required 
for self-similarity. As discussed at the end of the previous section. Figure 3 shows that for 
simulation C the power on A: > 80, i.e. on k larger than l/6th the Nyquist frequency of 
the PM grid, is artificially suppressed by 50% at the latest a with kni{a) = 7.25 (in units 
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of 2tt/L). The n = —2 simulation of Padmanabhan et al. (1995) is a PM simulation with 
240^ particles on a 720^ PM mesh, and used a staggered mesh scheme. Thus its spatial 
resolution is about twice superior to that of simulation C, but is still considerably worse 
than that of simulation B (since its PM grid spacing is > 2e for simulation B) which is 
shown in Figure 6 to lack the resolution to test self-similar scaling. 


4.1. Shape of the Nonlinear Spectrum 


Having shown that simulation A for the n = —2 spectrum follows self-similar scaling, 
we can characterize the spectrum at any time by the self-similar functional form of equation 
(§). The qualitative features of this shape can be seen in Figure 1 which shows the nonlinear 
and linear P(a, k) at hve different times, and in Figure 5 which shows the self-similar shape 
of the spectrum traced out by the measured spectra at eight different times. At sufficiently 
low {k/kni) the spectrum follows the linear shape oc k~‘^. At /c ~ kni the spectrum rises 
above the linear spectrum and then at still higher k it again approaches the shape oc k~‘^. 
The enhancement of the nonlinear spectrum relative to the linear spectrum is a factor of 
~ 1.5 at fc = kni, and approaches a factor ~ 5 at fc ^ 10 kni- It maintains this enhancement 
up to the highest k measured. 

The shape of P{k) measured at fc ^ 10 kni is in agreement with the prediction of the 
stable clustering ansatz. This ansatz relies on the assumption that the mean pair velocity 
in physical coordinates is zero. As shown by Davis & Peebles (1977) this leads to the 
prediction that the asymptotic form of P{k) at high-Zc is: P{a,k) oc a6/(5+n)^-6/(5+n)^ 
n = —2 this is: P oc a^k~‘^, in excellent agreement with the measured slope shown in Figure 
5. Note that once self-similarity is taken to be valid, the growth in a is fixed by the k~‘^ 
shape. Note also that for n = —2 the shape of the spectrum in the stable clustering regime 
is the same as that of the linear spectrum. 

The transition from the linear regime to the deeply nonlinear regime of stable clustering 
is poorly understood analytically, since perturbation theory breaks down at fc ~ kni- 
However a semi-empirical prescription for obtaining the nonlinear ^ for arbitrary initial 
spectra was proposed by Hamilton et al. (1991), and extended to the power spectrum 
by Peacock & Dodds (1994). Jain, Mo & White (1995) and Peacock & Dodds (1996) 
have refined the prescription to take into account a dependence on the initial spectrum. 
They provide fitting formulae for the dimensionless nonlinear spectrum at a given time, 
A(a, k) = 47r/c^P(a, k) in terms of the linear spectrum Ai(a, ki) = 47rfc^Piin(a, ki) evaluated 
at a smaller wavenumber ki as indicated. The two wavenumbers ki and k are related by 
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k = [1 + A{k)Y^^kL. One has to first obtain A{k) for given kL and using equation || 
below, and then substitute it into the relation above to get the value of k. The relation 
between A(a, k) and Ai(a, ki) is (Jain, Mo & White 1995): 


where the constant B{n) 


A{a,k) ^ ^ r Aj^(a, ki) 

B{n) [ B{n) \ ’ 

[(3 + n)/3]^'^, and the function <l> is 


<h(a;) = X 


'1 + 0.6x + x^- 0.2x^ - 

1 + 0.0037x3 ) 


(9) 


( 10 ) 


The formulae proposed by Peacock & Dodds (1996) give similar results as the equations 
above, though the nonlinear asymptote is somewhat lower for the n = —2 spectrum. 


The htting formula given above can be tested by using the nonlinear spectrum 
measured in our simulations. The results for the spectrum from simulation A are shown in 
Figure 7. The dimensionless spectrum A{a,k) is plotted against A 2 ,(a, fc^) for the output 
times shown in Figure 1. The solid curve shows the £t provided by equations (|^) and ([lOD . 
The dashed curve is obtained from the fitting formula of Peacock & Dodds (1996). The 
dotted lines show the slope of the linear and stable clustering predictions. The agreement 
between the N-body spectra and the solid curve is adequate, and is better than about 30% 
over almost the entire range plotted (note that the accuracy of the N-body spectrum at the 
smallest and largest amplitudes shown is about 20%). The N-body spectrum is enhanced 
relative to the fit at large amplitudes. This is in the same sense as the discrepancy noted 
by Jain, Mo & White (1995) who observed that the n = —2 spectrum has a slightly steeper 
nonlinear shape in the A{k) —A^lki) plot as compared to spectra with larger n. The results 
also show that the N-body spectra approach the stable clustering slope at amplitudes of 
A{k) > 50. 


5. Alternate Measures of Scaling 


The scaling properties of the matter distribution for scale-free conditions can be 
measured from any appropriate statistical measure. The power spectrum is the second 
moment of the density held, and is therefore one such statistic. The previous section 
demonstrated that it is only with very large simulation that the scaling of the power 
spectrum can be accurately tested. Higher order statistics are even harder to measure 
accurately. Of course every statistical measure which relates to the growth of perturbations 
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should measure the same scaling properties, but to measure them in N-body simulations 
for “difficult” spectra like n = —2, it is important to check more than one statistic. In 
this section we measure the scaling properties of the full density field by using a different 
approach. Rather that compute a statistic and then measure its scaling in time, we compute 
the fractional departure of the phase and amplitude from linear evolution mode by mode. 
The degree of nonlinearity thus measured is used to dehne characteristic nonlinear scales for 
the phase and amplitude. Thus by using both the phase and amplitude, and by using the 
nonlinearity of individual modes rather than an averaged quantity like the power spectrum, 
we are able to probe the scaling of the full density held more completely. 

We use the data for 5{k^ a) to compute the amplitude A and phase 0 given by, 

a) = A(fc, (11) 

A and 0 are the basic variables used for the scaling analysis in this section. In Figures 
8-10 the trajectories of the phase and amplitude of individual Fourier modes as a function 
of time are shown for the n = —2 simulation. Each of these hgures has four panels, and 
each panel shows the evolution of hve modes, chosen so that they represent a large range 
in k. These hgures give an idea of how individual modes evolve, in contrast to the regular 
behavior shown by the statistics computed from them. Even at relatively early times when 
most statistics obey linear behavior the amplitudes and phases can be seen to follow quite 
jagged paths, with the amplitude even showing negative growth for some time intervals. 
This suggests that there is more information to be mined than is provided by conventional 
statistics. 


5.1. Measuring the Phase Shift 


A technical problem in the measurement of the phase arises because conventionally the 
phases are dehned modulo (27r). The phase trajectories that result are shown in Figure 8, 
where all the phases he between —tt and tt. This is how the phases from N-body data have 
been computed in previous studies (Ryden & Gramann 1991; Scherrer, Melott & Shandarin 
1991; Suginohara & Suto 1991). However, following the phase trajectories in Figure 8 
makes it clear that this will have the effect of randomizing the phases at late times relative 
to their initial values, even if the actual growth were monotonic. This is because even a 
small change in the phase, 50, could cause it to be mapped to a value indicating a change 
of (27r — |50|). 



15 


If the trajectories could be obtained with arbitrarily small increments in a then such 
artihcial mappings could be un-done, and the phases plotted without constraining them 
between —n and n. Since outputs are available only at discrete values of a, there is a 
two-fold ambiguity in dehning the phase. Consider the phase values at two successive 
a’s for a given k: (pik^aj) and 0(A;, Oj+i), dehned in the usual way to lie between —tt 
and TT. Let dpj = 0(/c, Oj+i) — pik^aj). An alternate value for the phase at Oj+i is 
Oj+i) = (j){k,aj) ± (27r — |(50jD, where the sign in positive if Spj < 0 and vice versa. 
To choose between 0 and 0', we follow the trajectory of each mode, and at each snccessive 
valne of a, we dehne the phase by taking the magnitnde of the change in phase to be the 
smaller of \S(j)j\ and (27r — 1500). The resnlt is shown in Figure 9. As long as the typical 
changes in phase at snccessive times are less than tt, this procednre is a reasonable way of 
extending the range of 0. As we shall see, this considerably extends the degree of phase 
nonlinearity accessible to onr analysis. 


5.2. Scaling of the Phase and Amplitnde 


Following Ryden & Gramann (1991) and Gramann (1992), we dehne the following 
statistics as a measnre of the degree of nonlinear evolution. For the phase we dehne the 
mean deviation from the initial phase, 50(A;, a) = (|50(A;,a)|) = (|0(A;, a) — 0(A;,aj)|), where 
is the initial valne of a. The averages indicated are performed over the diherent modes 
within a shell in /c—space whose wavennmbers lie between [k — 0.5) and {k + 0.5). For the 
amplitude we simply measure the mean amplitnde (A(/c,a)) within each shell in /c—space. 
If we had used the phase trajectories as shown in Fignre 8 (i.e., dehned to he between — tt 
and tt), then at late times S(j){k,a) would have reached a maximum value 27r/3 — this 
corresponds to a distribntion of pik^a) that is nncorrelated with (f){k,ai). However, we hnd 
that 50 shows systematic growth well beyond 27r/3. Thus with the phase information that 
we have generated, previously unexplored aspects of phase evolution in the deeply nonlinear 
regime can be addressed. 


To analyze self-similar scaling, we dehne two characteristic wavennmbers. The hrst, 
denoted /Cc(a, 0c), is dehned by setting 50(/c, a) = 0c, where 0c is a constant. The second, 
denoted kc{a, Ac) is dehned using the amplitnde as follows: 


\A{k, a) — Ai{k, a)\ 
Ai{k,a) 


= Ac, 


( 12 ) 


where Ai(A;,a) is the linear solntion for A, and Ac is a dimensionless constant. Thus 
kc{a, Ac) is the wavenumber at which the fractional departnre of the amplitnde from the 
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linear solution is A,,. These statistics involve summing the magnitudes of the departures 
from linear behavior for each mode within a given fc—shell. Hence they probe the degree 
of nonlinearity more directly than if a statistic was computed hrst, and then its departure 
from the linear solution was calculated. 


The analytical prediction for the variance of the phase shift given in Paper I is 




1 2 


A-tt c 

= —a(Tfk'‘jdk,PUh 


(13) 


Thus the leading order solution for 5{k) involves a growing (and, for n < —1, divergent) 
phase shift, but there are no contributions to the amplitude at this order. With some 
further assumptions we were able to show analytically that the amplitude does not diverge 
provided n > —3, and should therefore show the standard self-similar scaling. We did not 
however obtain an analytical expression analogous to ([T^) for its growth. 


Figure 11 is a plot of log[/Cc(a, 4>c)] vs. a with n = —2, for 4 different values of (pc- Also 
shown in the plot are the scalings expected from self-similarity, k oc a“^, and the scaling 
resulting from the solution for 5(p given by equation (^), k oc a~^. This is the same as the 
velocity scaling referred to in Sections 2 and 3. The plots show that for 0^ = vr and 7r/2, the 
velocity scaling is closely followed, typically to better than 20%; but for the lower values 
(pc = 7 i/4: and tt/S, self-similar scaling is more closely followed at the same level of accuracy. 
We verihed that the trends did reflect a gradual transition in the scaling of kc{a, (pc) by 
plotting a larger range of (pc down to (pc = 7r/20. 


The dual scaling behavior shown by the phase shift can be interpreted as follows. For 
a given k, at early times as the phase just begins to depart from the linear solution (in 
which (p remains constant in time), its evolution is dominated by perturbative or other 
weakly nonlinear effects. These effects in general involve the coupling of a range of values 
of k', mostly in the vicinity of k' = A;, and obey the standard self-similar scaling. However 
at late times the phase shift is dominated by the bulk flow due to the longest waves in the 
box. The resulting phase change is given by equation ([T^) — it is a smooth function of a 
and /c, and therefore dominates the more stochastic, dynamical components of the phase 
change at late times. Thus the phase shows behavior that we can interpret as arising from 
a combination of the kinematical divergence, and a dynamical, nonlinear component which 
obeys the standard self-similar scaling. The former drives kc{a, (pc) to the kc oc a~^ scaling 
at late times, and the latter to the kc oc scaling at early times. 

Figure 13 shows kc{a, (pc) vs. a for the n = 0 spectrum. These results are from a P^M 
simulation performed by S. White, with 100^ particles and e = L/2500. The n = 0 spectrum 
shows only one behavior, the self-similar scaling, k oc This is expected as the linear 





17 


bulk velocity does not diverge, therefore the longest waves in the box do not dominate the 
phase shift at any time. 

For the amplitude scaling, Figure 12 shows a plot of log[A;c(a, Ac)] vs. a with n = —2, 
for Ac = 0.25, 0.5,1, 2. For sufficiently high k, all four curves closely follow the standard 
self-similar scaling, k oc a“^, to better than 20% accuracy. This is consistent with the 
results for the power spectrum, as one would expect since the power spectrum measures 
the variance of the amplitude. All the curves show a departure from the k oc scaling 
at low k, for k < 10, with the discrepancy exceeding a factor of 2 in the amplitude. This 
most likely indicates that the absence of power on modes with wavelengths larger than the 
box-size has slowed the growth of modes which would otherwise be enhanced by coupling to 
modes longer than the box. Thus the standard self-similar solution for n = —2 is obtained 
only on scales signihcantly smaller than the box-size. For n = 0 the self-similar scaling 
k oc is again shown convincingly in Figure 14. 

Our results are in partial disagreement with those of Gramann (1992). She found that 
for n = —1 in two dimensions (the analog of n = —2 in three dimensions), the standard 
self-similar scaling is broken for both the phase and amplitude. Our results for the phase 
scaling are consistent with hers, but the amplitude scalings are quite different: our results 
show good agreement with the scaling kc oc a“^, whereas hers agree with kc oc a~^. Since the 
statistics that we have measured are exactly the same as hers, it is difficult to explain the 
origin of the disagreement. It is conceivable that there are basic differences in the dynamics 
in two and three dimensions, but this is not reflected in the analytical results of Paper I. It 
is also possible that effects of the hnite box-size are more prominent in two dimensions. 


6. Conclusion 


In Section 2 we summarized the motivation for examining the self-similar scaling of 
scale free spectra, P{k) oc k^, for —3<n< —1. We have examined this issue through 
analytical techniques in Paper I and N-body techniques in this paper. In Paper I we 
found through two different approaches that the self-similar scaling of the amplitude of the 
density, and therefore of any measure of dynamical evolution, is preserved for —3 < n < 1. 

In this paper we have tested the scaling of n = —2 scale free simulations. In Section 
3 we discuss the issues of numerical resolution involved in such tests by comparing the 
resolution of three different simulations. We hnd that for the n = —2 spectrum we need the 
full resolution of our largest simulation, a P^M simulation with 256^ particles to measure 
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its scaling properties. Section 4 shows that self-similar scaling is verihed by the power 
spectrum measured from this simulation. Figures 4 and 5 demonstrate the existence of this 
scaling, and contrast it with an alternate scaling driven by the large-scale velocity held. We 
have compared our results with the recent work of Colombi, Bouchet & Hernquist (1995); 
Jain, Mo & White (1995); Lacey & Cole (1994); and Padmanabhan et ah (1995). We hnd 
that the shape of the nonlinear spectrum is well ht by the formulae proposed by Jain, Mo 
& White (1995) and Peacock & Dodds (1996), though the formulae slightly underestimate 
the spectrum in the strongly nonlinear regime. The shape of P{k) at high-/c agrees with the 
stable clustering prediction P oc k~^. 

The second test of the scaling properties of the density held made in Section 5 relies 
on the mode-by-mode evolution of the amplitude and phase. At sufficiently late times 
the phase shift obeys the velocity scaling, consistent with the solution found in Paper I. 
At early times the evolution of the phase shift is consistent with the standard self-similar 
scaling. Thus the phase shift arises from a combination of kinematical ehects due to large 
scale hows which dominate at late times, and genuine dynamical ehects which dominate 
at relatively early times — in the weakly nonlinear stage of its evolution. The scaling of 
the amplitude follows the standard self-similar form, except at wavenumbers /c ^ 10 (in 
units of 271/L) which we believe is a numerical limitation due to the hnite size of the box. 
This is consistent with the analytical results of Paper I and with the results for the power 
spectrum. 

In combination with the analytical analysis of Paper I, our results lead us to conclude 
with some conhdence that the self-similar evolution of the density contrast is preserved for 
—3 < n < —1. The kinematical interpretation for the scaling of the phase shift provides 
a useful guide to identifying statistics susceptible to such effects. The rms displacement 
of particle positions is an example of a statistic which would be dominated by the bulk 
motions from long wave modes for n < —1, and must therefore be used with caution as a 
measure of nonlinear evolution. 
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Fig. 1.— The power spectrum P{k) vs. k measured from simulation A at five different 
times. The solid curves show the N-body power spectrum for kni{a) = 58, 29,14.5, 7.25. The 
wavenumber k is in units of 27r/L. The dotted lines show the linear spectrum P oc a^k~‘^ 
at the same times. The N-body spectrum at low-k shows statistical fluctuations due to the 
small number of modes available to measure it. At high-Zc the spectrum is plotted up to 1/2 
the Nyquist frequency of the FFT grid, k = 512. 
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Fig. 2.— The effect of the FFT grid on P{k) at high-Zc. The two solid curves show the 
power spectrum computed from simulation B with different grid sizes: 512^ and 1024^. The 
spectrum measured with a 512^ grid is plotted with a dotted curve between k = 257 — 512, to 
show the artihcial features introduced by the grid used to interpolate the real space density. 
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Fig. 3.— The effect of limited resolution on the power spectrum. The solid curves show 
the power spectrum measured from simulation C, a PM simulation with the same number 
of particles as simulation A, a P^M simulation. The spectra are computed on a 512^ grid 
and are shown at the same output times as those of Figure 1. The dotted curves show the 
spectra from simulation A, computed using a 1024^ grid. Both spectra are plotted up to 1/2 
the Nyquist frequency of the FFT grid. A comparison of the two sets of spectra shows that 
as time evolves, the power in the PM simulation is suppressed on an increasingly large range 
of k. 
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Fig. 4.— The scaling of the power spectrum for the n = —2 simulation A. Simulation A is 
a P^M simulation with 256^ particles. The solid curves show the N-body power spectrum at 
eight different output times spanning kni{a) = 58 — 7.25. The dotted curves are all obtained 
by scaling the spectrum at the fourth from last output time: the wavenumber is scaled as 
k oc a(r)“", and the spectrum is scaled to maintain consistency with the linear spectrum. 
The left panel shows the results for the standard self-similar scaling with a = 2. The dotted 
curves agree very well with the solid curves except at the smallest and largest k owing to the 
hnite resolution of the simulation. The right panel shows the velocity scaling with a = 1. 
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Fig. 5.— The same power spectra as in Figure 4 are shown, but the scaling is tested a little 
differently. The power spectra and wavenumbers at the eight output times are scaled so that 
they would all he along a single curve if the scaling were correct. Again the standard self¬ 
similar scaling shown in the left panel works better than the alternate velocity scaling shown 
in the right panel. The lower dashed line shows the linear spectrum. The upper dashed line 
has arbitrary amplitude, but its slope is predicted by the stable clustering hypothesis. The 
N-body spectra asymptote to this slope at high-Zc. 
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Figure 6 


Fig. 6.— This is identical to Figure 5, but for simulation B which has 128^ particles and 
twice as large a spatial resolution as simulation A. The results shown in the two panels are 
not able to distinguish the different scalings. This hgure demonstrates the need for larger 
simulations, like the one shown in Figures 4 and 5, to effectively test the scaling properties 
of the n = —2 spectrum. 
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Fig. 7.— Comparison of N-body spectra with analytical fitting formulae. The dotted 
curves show the dimensionless power A(/c) = 47ik^P{k) plotted against the linear power 
on wavenumber as described in Section 4.1. The curves are obtained from the N-body 
spectrum of simulation A at the same times as in Figure 1. The solid curve is the htting 
formula of Jain, Mo & White (1995) as given in equation (PUD and the dashed curve is from 
the htting formula of Peacock & Dodds (1996). The upper dotted line shows the slope 
predicted by the stable clustering hypothesis, while the lower dotted line is the linear theory 
relation. 




Fig. 8.— Phase trajectories for the n = —2 simulation B. The trajectories of the phases 
(f>{k,T) of individual Fourier modes for n = —2 are plotted vs. a(r). The magnitudes k 
of the wavevectors are labeled in each panel; the full vectors were chosen as /c = (0,0,/c). 
Within each panel k increases in the following order: solid, dotted, dashed, long-dashed, 
dashed-dotted curves. The phase is dehned modulo 27r, and is therefore constrained to he 
between —tt and tt. At a = 0, the value of 0 at the earliest time has been plotted again to 
show the expected linear behavior. 


















-29- 



Fig. 9.— “Re-defined” phase trajectories for n = —2. The same trajectories as in Figure 8 
are plotted, but 0 has been re-de£ned so that it is no longer constrained between —tt and 
71 (notice the limits on the y-axis), as described in Section 5.1. In linear theory the phases 
do not change with time; significant departures from this can be seen in all but the lowest k 
modes. 
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Fig. 10.— Amplitude trajectories for n = —2. The trajectories of the amplitudes A{k,T) 
of individual Fourier modes for n = —2 are plotted vs. a, as in Figure 9 for the phase. 
Note that at early times A oc a: to check this all the curves have been joined to A = 0 at 
a = 0. Therefore the lowest value of a at which departures from a straight line occur shows 
nonlinear behavior. 
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Fig. 11.— Characteristic scales from the phase for n = —2. The scaling of characteristic 
wavenumber scales, kc vs. a is shown, as derived from the mean phase deviation. kc{a, 0c) 
is the value of k at which (|(50(/c, a) |) = (pc- The four solid curves correspond to the values 
of pc labeled at the top of the plot. Note that for high values of pc (the top 2 curves), the 
scaling closely agrees with kc cc shown by the upper dotted line. A transition towards 
the standard self-similar scaling kc oc a“^, shown by the lower dotted line, occurs for the two 
lower values of pc- 
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Fig. 12.— Characteristic scales from the amplitude for n = —2. The scaling of characteristic 
wavenumber scales derived from the departure of the amplitude from linear growth, kc{a, Ac) 
vs. a is shown for 4 different values of Ac (see equation ( p!2[) for the dehnition of kc{a, Ac)). 
All the curves have a slope close to the standard self-similar scaling, kc oc at high k. For 
k below about 10 the slope becomes shallower, probably due to the limitation of a hnite box. 
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Fig. 13.— Characteristic scales from the phase for n = 0. For n = 0 the scaling kc vs. a 
is shown using the mean phase deviation. Notice that in contrast to the n = —2 case, here 
only one behavior for all values of (be is evident: the self-similar scaling kc oc 
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Fig. 14.— Characteristic scales from the amplitude for n = 0. For n = 0 the scaling kc 
vs. a is shown using the mean amplitude. Again, as in Figure 13, the standard self-similar 
scaling is recovered. 




